#include "reduze.h"
#include <boost/graph/prim_minimum_spanning_tree.hpp>

scalar_products reduze_sp()
{
  using namespace cnst;
  scalar_products sp;
  sp.add(k1, i1, k1_i1);
  sp.add(k1, i2, k1_i2);
  sp.add(k1, i3, k1_i3);
  sp.add(k1, i4, k1_i4);
  sp.add(k1, i5, k1_i5);
  sp.add(k1, i6, k1_i6);

  sp.add(k2, i1, k2_i1);
  sp.add(k2, i2, k2_i2);
  sp.add(k2, i3, k2_i3);
  sp.add(k2, i4, k2_i4);
  sp.add(k2, i5, k2_i5);
  sp.add(k2, i6, k2_i6);

  sp.add(k3, i1, k3_i1);
  sp.add(k3, i2, k3_i2);
  sp.add(k3, i3, k3_i3);
  sp.add(k3, i4, k3_i4);
  sp.add(k3, i5, k3_i5);
  sp.add(k3, i6, k3_i6);

  sp.add(k4, i1, k4_i1);
  sp.add(k4, i2, k4_i2);
  sp.add(k4, i3, k4_i3);
  sp.add(k4, i4, k4_i4);
  sp.add(k4, i5, k4_i5);
  sp.add(k4, i6, k4_i6);

  sp.add(k5, i1, k5_i1);
  sp.add(k5, i2, k5_i2);
  sp.add(k5, i3, k5_i3);
  sp.add(k5, i4, k5_i4);
  sp.add(k5, i5, k5_i5);
  sp.add(k5, i6, k5_i6);

  sp.add(k6, i1, k6_i1);
  sp.add(k6, i2, k6_i2);
  sp.add(k6, i3, k6_i3);
  sp.add(k6, i4, k6_i4);
  sp.add(k6, i5, k6_i5);
  sp.add(k6, i6, k6_i6); 



  sp.add(k1, k1, k1_k1);
  sp.add(k1, k2, k1_k2);
  sp.add(k1, k3, k1_k3);
  sp.add(k1, k4, k1_k4);
  sp.add(k1, k5, k1_k5);
  sp.add(k1, k6, k1_k6);
  sp.add(k2, k2, k2_k2);
  sp.add(k2, k3, k2_k3);
  sp.add(k2, k4, k2_k4);
  sp.add(k2, k5, k2_k5);
  sp.add(k2, k6, k2_k6);
  sp.add(k3, k3, k3_k3);
  sp.add(k3, k4, k3_k4);
  sp.add(k3, k5, k3_k5);
  sp.add(k3, k6, k3_k6);
  sp.add(k4, k4, k4_k4);
  sp.add(k4, k5, k4_k5);
  sp.add(k4, k6, k4_k6);
  sp.add(k5, k5, k5_k5);
  sp.add(k5, k6, k5_k6);
  sp.add(k6, k6, k6_k6);

  sp.add(i1,i1, i1_i1);
  sp.add(i2,i2, i2_i2);
  sp.add(i3,i3, i3_i3);
  sp.add(i4,i4, i4_i4);
  sp.add(i5,i5, i5_i5);
  sp.add(i6,i6, i6_i6);

  sp.add(i1,i2, i1_i2);
  sp.add(i1,i3, i1_i3);
  sp.add(i1,i4, i1_i4);
  sp.add(i1,i5, i1_i5);
  sp.add(i1,i6, i1_i6);

  sp.add(i2,i3, i2_i3);
  sp.add(i2,i4, i2_i4);
  sp.add(i2,i5, i2_i5);
  sp.add(i2,i6, i2_i6);

  sp.add(i3,i4, i3_i4 );
  sp.add(i3,i5, i3_i5 );
  sp.add(i3,i6, i3_i6 );

  sp.add(i4,i5, i4_i5);
  sp.add(i4,i6, i4_i6);

  sp.add(i5,i6, i5_i6);

 

  //  sp.add(i1, i1, pow(m,2)); // A and B are orthogonal
  //  sp.add(k1, k1, 0); // A and C are orthogonal
  //  sp.add(k2, k2, 0); // A^2 = 4 (A has length 2)
  //  sp.add(k1, i1, k1_i1);
  // sp.add(k1,k2, k1_k2);
  //  sp.add(k2, i1, k2_i1);
  return sp;

}
/* changed to functor!!!
scalar_products sp_from_ex(ex in_ex)
{
  //  std::set<symbol> 
exset sym_set;
  for(const_preorder_iterator it = in_ex.preorder_begin();it!=in_ex.preorder_end();++it)
    {
      if(is_a<indexed>(*it))
        {
          ex ind_ex = *it;
          
          exmap repls;
          ind_ex.match(indexed(wild(),ind_ex.get_free_indices()),repls);
          cout<< "finding in "<< ind_ex<<" found "<<repls<<endl;
          ind_ex = wild().subs(repls);
          for(const_preorder_iterator s_it = ind_ex.preorder_begin();s_it!=ind_ex.preorder_end();++s_it)
            if(is_a<symbol>(*s_it))
              {
                cout<<"   "<<(*s_it)<<endl;
                sym_set.insert(ex_to<symbol>(*s_it));
              }
        }
    }
  cout<<sym_set<<endl;
  scalar_products sp;
  BOOST_FOREACH(ex sp_sym1,sym_set)
    {
      BOOST_FOREACH(ex sp_sym2,sym_set)
        {
          string str = ("(");
          sp.add(sp_sym1,sp_sym2,get_symbol(str+ex_to<symbol>(sp_sym1).get_name()+"."+ex_to<symbol>(sp_sym2).get_name()+")"));
        }
    }
  return sp;
}

*/
// find parts of scalar product by scalar product symbol
std::pair<ex,ex> sp_v(ex sp_ex)
{
  using namespace cnst;
  scalar_products rsp(reduze_sp());
  for(lst::const_iterator kit = km.begin();kit!=km.end();++kit)
    for(lst::const_iterator iit = im.begin();iit!=im.end();++iit)
      if(rsp.evaluate(*kit,*iit,wild())==sp_ex)return std::make_pair(*kit,*iit);
  for(lst::const_iterator kit = im.begin();kit!=im.end();++kit)
    for(lst::const_iterator iit = kit;iit!=im.end();++iit)
      if(rsp.evaluate(*kit,*iit,wild())==sp_ex)return std::make_pair(*kit,*iit);
  for(lst::const_iterator kit = km.begin();kit!=km.end();++kit)
    for(lst::const_iterator iit = kit;iit!=km.end();++iit)
      if(rsp.evaluate(*kit,*iit,wild())==sp_ex)return std::make_pair(*kit,*iit);
}

lst get_sp_lst()
{
  using namespace cnst;
  lst out_list;
  scalar_products rsp(reduze_sp());
  for(lst::const_iterator kit = km.begin();kit!=km.end();++kit)
    for(lst::const_iterator iit = im.begin();iit!=im.end();++iit)
      if(rsp.is_defined(*kit,*iit,wild()))out_list.append(rsp.evaluate(*kit,*iit,wild()));
  // for(lst::const_iterator kit = im.begin();kit!=im.end();++kit)
  // for(lst::const_iterator iit = kit;iit!=im.end();++iit)
  //   if(rsp.is_defined(*kit,*iit,wild()))out_list.append(rsp.evaluate(*kit,*iit,wild()));
  for(lst::const_iterator kit = km.begin();kit!=km.end();++kit)
    for(lst::const_iterator iit = kit;iit!=km.end();++iit)
      if(rsp.is_defined(*kit,*iit,wild()))out_list.append(rsp.evaluate(*kit,*iit,wild()));
  return out_list;
}

unsigned int ID(std::vector<int> vec)
{
  int id_out = 0;
  int ji = 0;
  BOOST_FOREACH(int nu,vec)
    {
      if(nu>0) id_out+=ex_to<GiNaC::numeric>(pow(2,ji)).to_int();
      ji++;
    }
  return id_out;
}
// summ of numerator powers
unsigned int s_nu(std::vector<int> vec)
{
  int id_out = 0;
  BOOST_FOREACH(int nu,vec)
    {
      if(nu<0) id_out-=nu;
    }
  return id_out;
}
unsigned int r_nu(std::vector<int> vec)
{
  int id_out = 0;
  BOOST_FOREACH(int nu,vec)
    {
      if(nu>0) id_out+=nu;
    }
  return id_out;
}
unsigned int t_nu(std::vector<int> vec)
{
  int id_out = 0;
  BOOST_FOREACH(int nu,vec)
    {
      if(nu>0) id_out++;
    }
  return id_out;
}

std::string nu_str(std::vector<int> vec)
{
  std::stringstream str_out;
  str_out<<"{";
  for(std::vector<int>::iterator it = vec.begin();it!=vec.end();++it)
    {
      if(vec.end() == boost::next(it))
        str_out<<*it;
      else
        str_out<<*it<<",";
    }
  str_out<<"}";
  return str_out.str();
}
// function splitting integral on integrals with different set of propagators
lst split_int(ex in_ex)
{
    //bubblesii::prop_pow_helper_t my1;
  lst out_lst;
  ex  int_set = in_ex.expand().collect(props(in_ex),true);
  cout<<"int_set "<<int_set<<endl;
  //  massert("huynya!!!! ",is_a<add>(int_set)||is_a<mul>(int_set));
   BOOST_ASSERT_MSG(is_a<add>(int_set)||is_a<mul>(int_set), "Integral is not a MUL and not a SUM");
  if(is_a<add>(int_set))
    for(const_iterator it = int_set.begin();it!=int_set.end();++it)
      {
        //cout<<"pr_po "<<my1(*it)<<endl;
        out_lst.append(*it);
        //      get_no_prop_part(*it);
        //cout<<"propro "<< props(*it)<<endl<<endl;
      }
  else if (is_a<mul>(int_set))
    out_lst.append(int_set);
  //  else throw std::logic_error("Integral is not a MUL and not a SUM");
  return out_lst;
}

// function splitting integral on integrals with different set of propagators in terms of P_i in mincer
lst split_int_P(ex int_set)
{
    //bubblesii::prop_pow_helper_t my1;
  lst out_lst;
  //ex  int_set = in_ex.expand().collect(props(in_ex),true);
  //cout<<"int_set "<<int_set<<endl;
  //  massert("huynya!!!! ",is_a<add>(int_set)||is_a<mul>(int_set));
   BOOST_ASSERT_MSG(is_a<add>(int_set)||is_a<mul>(int_set), "Integral is not a MUL and not a SUM");
  if(is_a<add>(int_set))
    for(const_iterator it = int_set.begin();it!=int_set.end();++it)
      {
        //cout<<"pr_po "<<my1(*it)<<endl;
        out_lst.append(*it);
        //      get_no_prop_part(*it);
        //cout<<"propro "<< props(*it)<<endl<<endl;
      }
  else if (is_a<mul>(int_set))
    out_lst.append(int_set);
  //  else throw std::logic_error("Integral is not a MUL and not a SUM");
  return out_lst;
}


// !!! ONLY MINCER function, M=0!!!!
// only reducible numerators!!!
exmap  sp_in_prop(lst sp_lst_in,lst p_lst)
{
  try{
    symbol x0("x0");
    ex PX = x0;
    int cntr = 1;
    symbol idx_dummy("idx");
    string str;
    ex m_p1,m_p2;
    ex p_p1,p_p2;
    ex j_p;
    ex PX_symb = x0;
    lst x_lst;
  
    // testin sp in current prop
    lst sp_lst(get_sp_lst());
    cout<<"SPLST_GET"<<sp_lst<<endl;

    x_lst.append(x0);
    lst sp_lst_list;
    for(lst::const_iterator it1 = p_lst.begin();it1!=p_lst.end();++it1)
      for(lst::const_iterator it2 = it1;it2!=p_lst.end();++it2)
        {

          str = "x" + lexical_cast<string>(cntr);
          symbol x_sym(str);
          cout<<"Prop = 1/("<<ex_to<bubblesii::prop>(*it1).op(2)<<endl;
          m_p1 =  ex_to<bubblesii::prop>(*it1).op(0);
          m_p2 =  ex_to<bubblesii::prop>(*it2).op(0);
          j_p = ex_to<bubblesii::prop>(*it1).op(1);
          p_p1 = ex_to<bubblesii::prop>(*it1).op(2);
          p_p2 = ex_to<bubblesii::prop>(*it2).op(2);
          cout<<m_p1<<"  "<<p_p1<<"  "<<p_p2<<endl;

          lst N_has_sp_lst;
          ex in_ex_ex = (x_sym*(indexed(p_p1,varidx(idx_dummy,4)) *  indexed(p_p2,varidx(idx_dummy,4).toggle_variance()))).expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
      
          for(lst::const_iterator it_sp = sp_lst.begin();it_sp!=sp_lst.end();++it_sp)
            {
              if(in_ex_ex.has(*it_sp))N_has_sp_lst.append(*it_sp);
            }
          cout<<"Listochek "<<N_has_sp_lst<<endl;
          if(true)//!sp_lst_list.has(N_has_sp_lst.sort()))
            {
              x_lst.append(x_sym);
              sp_lst_list.append(N_has_sp_lst.sort());
              PX+=(x_sym*(indexed(p_p1,varidx(idx_dummy,4)) *  indexed(p_p2,varidx(idx_dummy,4).toggle_variance()))).expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
              PX_symb+=(*it1)*x_sym*(*it2); /// why division?????
              //D_symb*=p_sym;
              cntr++;

            }
          else cout<<"Listo not uni"<<endl;
        }
    cout<<"Listochek end"<<sp_lst_list.unique()<<endl;
    // END OF PX - linear combination of all scalar products


    cout<<"PPXX  "<<PX_symb<<endl;
    exmap sp_P_subs_map;
    lst PX_has_sp_lst;
    // filling list of sp in denominator
    for(lst::const_iterator it = sp_lst.begin();it!=sp_lst.end();++it)
      if(PX.has(*it)) PX_has_sp_lst.append(*it);
    cout<<"PXHSP  only in den: "<<PX_has_sp_lst<<endl;
    for(lst::const_iterator nit = sp_lst_in.begin();nit!=sp_lst_in.end();++nit)
      {
        // ex (*nit) = N_has_sp_lst.op(0);
        ex PX_sub_part = 0;
        lst eqns;
        for(lst::const_iterator it = PX_has_sp_lst.begin();it!=PX_has_sp_lst.end();++it)
          {
            if((*nit) == *it) 
              {
                eqns.append(PX.coeff(*it,1)==1);
                PX_sub_part += PX.coeff(*it,1)*(*it);           
              }
            else
              {
                eqns.append(PX.coeff(*it,1)==0);
                PX_sub_part += PX.coeff(*it,1)*(*it);           
              }
          }
        eqns.append((PX-PX_sub_part).collect(sp_lst) ==0);


        cout<<"EQNS_d var: "<<*nit<<"   "<<eqns<<" X:"<<x_lst<<
          /*"P_map "<<P_map<<"    PX      "<<PX<<"   sp__has  "<<irreducible_sp_lst<<"  scalar part  "<<(PX-PX_sub_part).collect(sp_lst)<<"  eqns   "<<eqns<< "    result    "<<
           */lsolve(eqns,x_lst)<<endl;                    
        sp_P_subs_map[*nit] = PX_symb.subs(lsolve(eqns,x_lst));
      }
    return sp_P_subs_map;
  }
  catch(std::exception &p)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<p.what()<<endl;
      cout<<"*********************************************"<<endl;
    }
}

// return list of pairs <coeff,integral>, where integral has only propagators
std::list<expair>  N_to_prop(ex in_ex)
{
  std::list<expair> out_lst;
  in_ex = in_ex.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  cout<<"spi "<<props(in_ex)<<endl<<endl;
    lst sp_lst(get_sp_lst());
    lst N_has_sp_lst;
    exmap sp_P_subs_map;
    for(lst::const_iterator it = sp_lst.begin();it!=sp_lst.end();++it)
      {
        if(in_ex.has(*it))N_has_sp_lst.append(*it);
    }
    if(N_has_sp_lst.nops())
      {
        exmap prop_over_sp(sp_in_prop(N_has_sp_lst,props(in_ex)));
        cout<<"spi "<<N_has_sp_lst<<endl<<endl;
        lst ilst =  split_int_P(in_ex);
        //bubblesii::prop_pow_helper_t int_wo_prop_pow;
        cout<<ilst<<endl;
        for(lst::const_iterator it = ilst.begin();it!=ilst.end();++it)
          {
            //        cout<<"spi "<<((it->subs(prop_over_sp)).expand()).collect(props(*it),true)<<endl;
            cout<<"spi  C= "<<get_no_prop_part(*it)<<"   I =  "<</*int_wo_prop_pow*/((*it)/get_no_prop_part(*it))<<endl<<endl;
            out_lst.push_back(expair(get_no_prop_part(*it),/*int_wo_prop_pow*/((*it)/get_no_prop_part(*it))));
          }
      }
    else
      {
        //exmap prop_over_sp(sp_in_prop(N_has_sp_lst,props(in_ex)));
        // cout<<"spi "<<N_has_sp_lst<<endl<<endl;
        cout<<"esplit :"<<in_ex.expand().collect(props(in_ex),true)<<endl;
        lst ilst =  split_int(in_ex.expand().collect(props(in_ex),true));
        cout<<ilst<<endl;
        //bubblesii::prop_pow_helper_t int_wo_prop_pow;
        for(lst::const_iterator it = ilst.begin();it!=ilst.end();++it)
          {
            //        cout<<"spi "<<((it->subs(prop_over_sp)).expand()).collect(props(*it),true)<<endl;
            cout<<"spi  C= "<<get_no_prop_part(*it)<<"   I =  "<</*int_wo_prop_pow*/((*it)/get_no_prop_part(*it))<<endl<<endl;
            out_lst.push_back(expair(get_no_prop_part(*it),/*int_wo_prop_pow*/((*it)/get_no_prop_part(*it))));
          }
      }
    return out_lst;
}

exmap subs_1_map(lst in_lst)
{
  exmap out_map;
  for(lst::const_iterator it = in_lst.begin();it!=in_lst.end();++it)
    out_map[*it] = 1;
  return out_map;
}

// return list of pairs <coeff,integral>, where integral has only propagators
// used P_map map to subs P momentums
std::list<expair>  N_to_mincer_P(ex in_ex,exmap P_map,exmap prop_map)
{
  std::list<expair> out_lst;
  in_ex = in_ex.expand();
  ex P_ex = 0;
  // subs k->P_i from mincer
 in_ex = in_ex.subs(P_map);
 // instance of sp symbol generator
  sp_from_ex_functor sp_from_ex;
  // convert bubblessi props to indexed expressions
  //  bubblesii_prop_to_idx_prop do_bub(sp_from_ex);
  //in_ex  = do_bub(in_ex);
 
  cout<<" sub_bef "<<in_ex<<endl;


  //  cout<<" sub_bef bub "<<do_bub(in_ex)<<endl;
//  in_ex = in_ex.normal();
   ex px_tmp2 = in_ex;
   cout<<" sub_aft "<<in_ex.expand(expand_options::expand_indexed).simplify_indexed(sp_from_ex(px_tmp2))<<endl;
   in_ex = in_ex.expand(expand_options::expand_indexed).simplify_indexed(sp_from_ex(px_tmp2));

   lst collect_lst(sp_from_ex.get_sp_symbols_lst());
   for(exmap::const_iterator it = prop_map.begin();it!=prop_map.end();++it)
     {
       collect_lst.append(1/(it->second));
     }
   cout<<" collect_lst: "<<collect_lst<<endl;
      cout<<"sub_aft_set " <<in_ex.expand().collect(collect_lst,true)<<endl;
      in_ex = in_ex.expand().collect(collect_lst,true);

   //   cout<<"sub_aft_set " <<in_ex.expand().collect(cnst::me,true)<<endl;
   //      in_ex  = in_ex.expand(expand_options::expand_indexed).simplify_indexed(sp_from_ex(px_tmp));

 
     lst sp_lst(get_sp_lst());
    lst N_has_sp_lst;
    exmap sp_P_subs_map;
    
    
      {
              
        lst ilst =  split_int(    in_ex);
        //bubblesii::prop_pow_helper_t int_wo_prop_pow;
        cout<<ilst<<endl;
        for(lst::const_iterator it = ilst.begin();it!=ilst.end();++it)
          {
            //        cout<<"spi "<<((it->subs(prop_over_sp)).expand()).collect(props(*it),true)<<endl;
            cout<<"spi  C= "<<it->subs(subs_1_map(collect_lst))<<"   I =  "<<((*it)/(it->subs(subs_1_map(collect_lst))))<<endl<<endl;
            out_lst.push_back(expair(it->subs(subs_1_map(collect_lst)),((*it)/(it->subs(subs_1_map(collect_lst))))));
          }
      }
    
    return out_lst;
}


// map between input graph edges and MINCER symbols for momentums
std::pair<exmap,bool> find_topo(Graph& G)
{
  try{
    exmap out_map;
  Graph_topo topo_vec (get_topo_v_e(num_vertices(G),num_edges(G)));
  cout<<"Find topo started"<<endl;
  property_map<Graph, vertex_index_t>::type 
    v1_index_map = get(vertex_index, G);
  std::vector<graph_traits<Graph>::vertex_descriptor> f(num_vertices(G));
  Graph::edge_descriptor edis;
  Graph::edge_iterator ei,ei_end;
  Graph::out_edge_iterator oei,oei_end;
  bool bexist;


  BOOST_FOREACH(Graph_topo_pair graph_topo,topo_vec)
    {
      if( isomorphism
          (G, graph_topo.first, isomorphism_map
           (make_iterator_property_map(f.begin(), v1_index_map, f[0]))))
        {

          for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
            G[*ei].visited = false;
          for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
            {
              //if(has_par_edges(*ei,G))
                {
                  if(!G[*ei].visited)
                      {
                        cout<<"ADD_PAIR EDGE(G): "<<*ei<<" v:s: "<<source(*ei,G)<<",t: "<<target(*ei,G)<<"iso: ( "<<f[source(*ei,G)]<<","<<f[target(*ei,G)]<<")"<<endl;
                        std::vector<Graph::edge_descriptor> vec_G = edges_u_v(source(*ei,G),target(*ei,G),G);
                        std::vector<Graph::edge_descriptor> vec_iso = edges_u_v(f[source(*ei,G)],f[target(*ei,G)],graph_topo.first);
                        cout<<"iso g "<<endl;
                        BOOST_FOREACH(Graph::edge_descriptor eee,vec_iso)
                          {
                            cout<<eee<<endl;
                          }

                        if(vec_G.size() == vec_iso.size())
                          {
                            for(int i = 0;i<vec_G.size();i++)
                              {
                                out_map[props
                                        (G[vec_G[i]].gexpr).op(0)] = graph_topo.first[vec_iso[i]].gexpr;
                                cout<<"map: "<<G[vec_G[i]].gexpr<<"  , "<<graph_topo.first[vec_iso[i]].gexpr<<endl;
                                G[vec_G[i]].visited = true;
                              }
                          }
                        else throw std::logic_error("ERROR: different number of parallel edges between isomorph vertices");
                      }

                }
                /* else
                {
                  tie(edis,bexist)=edge(f[source(*ei,G)],f[target(*ei,G)],graph_topo.first);
                  if(bexist) out_map[*ei] = graph_topo.first[edis].gexpr;
                  else throw std::logic_error("ERROR: No such edge in isomorphism!!!");
                  }*/
            }
          std::cout << "isomorph to  topology "<<graph_topo.second<<std::endl;
          return std::make_pair(out_map,true);
        }
      
    }
  return std::make_pair(out_map,false);


}
  catch(string& what)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      cout<<"*********************************************"<<endl;
    }
}


// map between input graph edges momentums  and MINCER symbols for momentums
// prop_map,mom_map,found
tuple<exmap,exmap,string,bool> find_topo_momentum(Graph& G)
{
  try{
    exmap out_map;
    exmap out_mom_map;
  Graph_topo topo_vec (get_topo_v_e(num_vertices(G),num_edges(G)));
  cout<<"Find topo started"<<endl;
  property_map<Graph, vertex_index_t>::type 
    v1_index_map = get(vertex_index, G);
  std::vector<graph_traits<Graph>::vertex_descriptor> f(num_vertices(G));
  Graph::edge_descriptor edis;
  Graph::edge_iterator ei,ei_end;
  Graph::out_edge_iterator oei,oei_end;
  bool bexist;


  BOOST_FOREACH(Graph_topo_pair graph_topo,topo_vec)
    {
      if( isomorphism
          (G, graph_topo.first, isomorphism_map
           (make_iterator_property_map(f.begin(), v1_index_map, f[0]))))
        {

          for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
            G[*ei].visited = false;
          for(tie(ei,ei_end) = edges(G);ei!=ei_end;++ei)
            {
              //if(has_par_edges(*ei,G))
                {
                  if(!G[*ei].visited)
                      {
                        cout<<"ADD_PAIR EDGE(G): "<<*ei<<" v:s: "<<source(*ei,G)<<",t: "<<target(*ei,G)<<"iso: ( "<<f[source(*ei,G)]<<","<<f[target(*ei,G)]<<")"<<endl;
                        std::vector<Graph::edge_descriptor> vec_G = edges_u_v(source(*ei,G),target(*ei,G),G);
                        std::vector<Graph::edge_descriptor> vec_iso = edges_u_v(f[source(*ei,G)],f[target(*ei,G)],graph_topo.first);
                        cout<<"iso g "<<endl;
                        BOOST_FOREACH(Graph::edge_descriptor eee,vec_iso)
                          {
                            cout<<eee<<endl;
                          }
                        BOOST_ASSERT_MSG(vec_G.size() == vec_iso.size(),"different number of vertices in isomorph graphs");
                        cout<<"goot "<<vec_G.size()<<endl;

                        print(G);
                        if(vec_G.size() == vec_iso.size())
                          {
                            for(int i = 0;i<vec_G.size();i++)
                              {
                                cout<<"goot "<<G[vec_G[i]].gexpr<<endl;
    cout<<"goot "<< graph_topo.first[vec_iso[i]].gmomentum<<endl;
                                out_map[props(G[vec_G[i]].gexpr).op(0)] = graph_topo.first[vec_iso[i]].gexpr;
                                out_mom_map[ex_to<bubblesii::prop>(props(G[vec_G[i]].gexpr).op(0)).op(2)] = graph_topo.first[vec_iso[i]].gmomentum;
                    
                                cout<<"map: "<<G[vec_G[i]].gexpr<<"  , "<<graph_topo.first[vec_iso[i]].gexpr<<endl;
                                G[vec_G[i]].visited = true;
                              }
                          }
                        else throw std::logic_error("ERROR: different number of parallel edges between isomorph vertices");
                      }

                }
                /* else
                {
                  tie(edis,bexist)=edge(f[source(*ei,G)],f[target(*ei,G)],graph_topo.first);
                  if(bexist) out_map[*ei] = graph_topo.first[edis].gexpr;
                  else throw std::logic_error("ERROR: No such edge in isomorphism!!!");
                  }*/
            }
          std::cout << "isomorph to  topology "<<topo_map[graph_topo.second]<<std::endl;
          return make_tuple(out_map,out_mom_map,string(topo_map[graph_topo.second]),true);
        }
      
    }
  return make_tuple(out_map,out_mom_map,string("0"),false);


}
  catch(std::exception& p)
    {
      //      cout<<"*********************************************"<<endl;
      // cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      //cout<<"*********************************************"<<endl;
      throw std::logic_error(std::string("In function \"find_topo_momentum\":\n |___> ")+p.what());
    }
}




void mincer(Graph G,dia_ext_legs legs)
{
  try{

    /*  !!!prim tree!!!
std::vector < graph_traits < Graph >::vertex_descriptor >
    p(num_vertices(G));

  prim_minimum_spanning_tree(G, &p[0]);


  for (std::size_t i = 0; i != p.size(); ++i)
    if (p[i] != i)
      std::cout << "parent[" << i << "] = " << p[i] << std::endl;
    else
      std::cout << "parent[" << i << "] = no parent" << std::endl;
    */
  

 
    exmap pro_map;
    exmap mom_map;
    bool topo_b_f;
    Graph g_in(G);
    exmap subs_i1_eq_Q;
    subs_i1_eq_Q[cnst::i1] = get_symbol("Q");
    graph_ex_subs(subs_i1_eq_Q,g_in);
    print(g_in);
    string topo_str;
    tie(pro_map,mom_map,topo_str,topo_b_f) =  find_topo_momentum(g_in);
    cout<<"suuuubs "<<map_to_subs(mom_map)<<endl;
    //    assert(false);
  ex tmpex;
  if(topo_b_f)
    {
      for(exmap::iterator it = mom_map.begin();it!=mom_map.end();++it)
        {
          cout<<it->first<<"   ------   "<<it->second<<endl;
          //tmpex =  props(it->first).op(0);

          //it->first = tmpex;
       }
     
      cout<<"MMMAAAPPP "<<mom_map<<endl;
    }
  else throw std::logic_error("ERROR: Topology not found");



  //  pro_map[cnst::i1_i1] = get_symbol("(Q.Q)");

  //topology_end
  std::stringstream out_frm;
  out_frm <<"out/dia"<<get_property(g_in,graph_num)<<".frm";
  std::ofstream form_out(out_frm.str().c_str());
 form_out<<"#-\n#define SCHEME \"0\"\n#include mincer2.h\n#include afdc.h\n.global\n";

 // call root_projector
Graph fix_G =  fix_rl(g_in);

 ex in_norm_expr = ex_from_graph(fix_G);//.simplify_indexed();

 cout<<"tracewith "<<in_norm_expr<<endl<<endl;

 exvector gmunu = in_norm_expr.get_free_indices();
 //******* vacuum polarization way
 BOOST_ASSERT_MSG(gmunu.size()==2,"Working on vacuum polarization => only 2 free indices");
 cout<<gmunu<<delta_tensor(ex_to<varidx>(gmunu[0]).toggle_variance(),ex_to<varidx>(gmunu[1]).toggle_variance())<<endl;

 in_norm_expr = (in_norm_expr*lorentz_g(ex_to<varidx>(gmunu[0]).toggle_variance(),ex_to<varidx>(gmunu[1]).toggle_variance())/3).simplify_indexed();//.coeff(cnst::me,0);
 //*******************************
 in_norm_expr = trace_ng(in_norm_expr);
 cout<<"trac trace "<< in_norm_expr<<endl;

 
  cout<<"traceless "<<" <--rls "<<in_norm_expr<<endl;
  //  std::list<expair> c_i_lst = N_to_prop(in_norm_expr);
  std::list<expair> c_i_lst = N_to_mincer_P(in_norm_expr.subs(pro_map),map_to_subs(mom_map),pro_map);


  cout<<"size "<<c_i_lst<<endl;
  int cn = 0;
  BOOST_FOREACH(expair c_i,c_i_lst)
    {
      form_out<<"Global indiaC"<<cn<<"="<<c_i.rest.subs(pro_map)<<";\n"
              <<".store\n";
      form_out<<"Global indiaD"<<cn<<"="<<c_i.coeff.subs(pro_map)<<";\n"
              <<"Multiply ep^3;\n"
              <<"#call integral("<<topo_str<<")\n"
        //              <<"#include id"<<get_property(G,graph_num)<<"mapping.h\n"
              <<".store\n";
      form_out<<"Global indiaE"<<cn<<"="<<"indiaC"<<cn<<"*"<<"indiaD"<<cn<<";\n"
              <<"Print;\n"
              <<".store\n";
      cn++;
    }

  form_out<<"Local oex = ";
  for(int tc = cn-1;tc>0;tc--)
    {
      form_out<<"indiaE"<<tc<<"+ ";
    }
  form_out<<"indiaE0;\n"
          <<"Bracket ep;"
          <<"Print +s;\n"
          <<".end\n";


  }
  catch(std::exception &p)
    {
      // cout<<"******************************************************"<<endl;
      //cout<<">>>  ERROR:  "<<setw(40)<<std::internal<<p.what()<<endl;
      //cout<<"******************************************************"<<endl;
      throw std::logic_error(std::string("In function \"MINCER\":\n |___> ")+p.what());
     }

}

void reduze(Graph& G)
{
try
  {
  //  cout<<"reduze expr to work with"<<ex_from_graph(G).normal().expand(expand_options::expand_indexed).simplify_indexed(reduze_sp())<<"  sp  "<<sp_v(cnst::k1_k3).first<<" "<<sp_v(cnst::k1_k3).second<<endl;
  ex in_norm_expr = ex_from_graph(G).normal();

  // Topo.in filling
  std::ofstream topo_in("Topo.in_");
  topo_in<<"LoopMomenta";
  for(lst::const_iterator it = cnst::km.begin();it!=cnst::km.end();++it)
    if(in_norm_expr.has(*it))topo_in<<" "<<*it;
  topo_in<<";"<<endl<<endl;
  topo_in<<"ExtMomenta";
  for(lst::const_iterator it = cnst::im.begin();it!=cnst::im.end();++it)
    if(in_norm_expr.has(*it))topo_in<<" "<<*it;
  topo_in<<";"<<endl<<endl
         <<"Dimension d;"<<endl<<endl;

  std::stringstream symbols_str,propagators_str;
  symbols_str<<"Symbol";
  ex  N = in_norm_expr.numer().expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  ex D =  in_norm_expr.denom();//.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  cout<<"reduze num  "<<N<<"   denom  "<<D<<endl;  
  Graph::edge_iterator ei,ei_end;
  exmap repls;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    if(D.match(wild(2)*(indexed(G[*ei].gmomentum,varidx(wild(1),4))*indexed(G[*ei].gmomentum,varidx(wild(1),4,true)) -wild(3)),repls )) cout<< "has momentum   "<<G[*ei].gmomentum<<"   "<<repls<<endl;
    else  if(D.match(wild(2)*(indexed(wild(6),varidx(wild(1),4))*indexed(wild(6),varidx(wild(1),4).toggle_variance()) +wild(3)),repls )) cout<< "has (-) momentum   "<<G[*ei].gmomentum<<"   "<<repls<<endl;
    else cout<< "no momentum   "<<G[*ei].gmomentum<<endl;
  symbol P0("P0");
  symbol x0("x0");
  lst P_lst;
  //P_lst.append(P0);
  exmap P_map;
  //P_map[P0] = 1;
  lst x_lst;
  x_lst.append(x0);
  ex factor = 1;
  int cntr = 0;
  string str;
  ex  PX;   // Denominator in form P_i*X_i, where P_i = (q_i^2-m_i^2)
  PX = x0;
  ex PX_symb = x0; // where P_i are symbols
  ex D_symb = 1;
  exset mass_set; //set for symbols storage, no duplicated symbols
  for(const_iterator it = D.begin();it!=D.end();++it)
    {
      if( is_a<GiNaC::numeric>(*it))
        factor *=(*it);
      else
        {
          ex prop_ex;
          exmap prop_repls;
          ex found_mass(wild(1));
          if(it->match(wild()-pow(wild(1),2),prop_repls))
            {
              found_mass = found_mass.subs(prop_repls);
              cout<<prop_repls<<" momentum "<<*it<<endl;
              prop_ex = *it;
              cntr++;
            }
          else if(it->match(wild()+pow(wild(1),2),prop_repls))
            {
              found_mass = found_mass.subs(prop_repls);
              cout<<prop_repls<<" (-) momentum "<<*it<<endl;
              factor*=(-1);
              prop_ex = (*it)*(-1);
              cntr++;
            }
          else throw  string("Bad devision in rational expression");
          //          cout<<prop_ex<<"  wildcard "<<indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true))-pow(wild(2),2)<<endl;
          // extracting mass and momentum from propagators
          mass_set.insert(found_mass);
          exset found_set;
          exmap find_match;
          if(prop_ex.find(indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true)),found_set))
            {
            if(found_set.size() == 1)
              {
                if( found_set.begin()->match(indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true)),find_match))
                  propagators_str<<"Propagator  "<<ex(wild()).subs(find_match)<<" "<<found_mass<<";"<<endl;
                  //  cout<<"momrepls "<<*(found_set.begin())<<"  Propagator  "<<ex(wild()).subs(find_match)<<" "<<found_mass<<";"<<endl;
              }
            else throw string("Bad denominator expression");
            }
          str = "P" + lexical_cast<string>(cntr);
          symbol p_sym(str);
          P_lst.append(p_sym);
          P_map[p_sym] = prop_ex.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
          str = "x" + lexical_cast<string>(cntr);
          symbol x_sym(str);
          x_lst.append(x_sym);
          PX+=x_sym*P_map[p_sym];
          PX_symb+=x_sym*p_sym;
          D_symb*=p_sym;

        }
    }
  topo_in<<"Symbol";
  for(exset::iterator it = mass_set.begin();it!=mass_set.end();++it)
    topo_in<<" "<<*it;
  
  topo_in<<";"<<endl<<endl
         <<propagators_str.str()<<endl;
  N *=factor;

  lst sp_lst(get_sp_lst());
  lst N_has_sp_lst;
  lst PX_has_sp_lst;
  lst irreducible_sp_lst;
  exmap sp_P_subs_map;
  for(lst::const_iterator it = sp_lst.begin();it!=sp_lst.end();++it)
    {
      if(N.has(*it))N_has_sp_lst.append(*it);
      if(PX.has(*it))PX_has_sp_lst.append(*it);
      if(N.has(*it)&&(!PX.has(*it)))irreducible_sp_lst.append(*it);
    }
  for(lst::const_iterator nit = N_has_sp_lst.begin();nit!=N_has_sp_lst.end();++nit)
    {
      // ex (*nit) = N_has_sp_lst.op(0);
      ex PX_sub_part = 0;
      lst eqns;
      for(lst::const_iterator it = PX_has_sp_lst.begin();it!=PX_has_sp_lst.end();++it)
        {
          if((*nit) == *it) 
            {
              eqns.append(PX.coeff(*it,1)==1);
              PX_sub_part += PX.coeff(*it,1)*(*it);           
            }
          else
            {
              eqns.append(PX.coeff(*it,1)==0);
              PX_sub_part += PX.coeff(*it,1)*(*it);           
            }
        }
      eqns.append((PX-PX_sub_part).collect(sp_lst) ==0);


      //  cout<<"P_map "<<P_map<<"    PX      "<<PX<<"   sp__has  "<<irreducible_sp_lst<<"  scalar part  "<<(PX-PX_sub_part).collect(sp_lst)<<"  eqns   "<<eqns<< "    result    "<<lsolve(eqns,x_lst)<<endl;                    
      sp_P_subs_map[*nit] = PX_symb.subs(lsolve(eqns,x_lst));
    }
  
  ex the_int = N.subs(sp_P_subs_map);
  for(lst::const_iterator it = P_lst.begin();it!=P_lst.end();++it)
    the_int/=(*it);
  ex distrib_int = the_int.expand().collect(P_lst,true);


  std::list<std::vector<int> > nu_vec_lst;
  std::set<unsigned int> ID_set;
  unsigned int r_max = 0;
  unsigned int s_max = 0;
  // output file for integrals
  std::ofstream MyIntegrals("afdc_Integrals");
  for(const_iterator it = distrib_int.begin();it!=distrib_int.end();++it)
    {
      std::vector<int> nu_vec;
      ex tmp_ex2 = *it;
      for(lst::const_iterator pit = P_lst.begin();pit!=P_lst.end();++pit)
          {
                nu_vec.push_back(-tmp_ex2.degree(*pit));
                //cout<<-tmp_ex2.degree(*pit);
                //MyIntegrals<<-tmp_ex2.degree(*pit);
                //               cout<<-tmp_ex2.degree(*pit)<<", ";
                // MyIntegrals<<-tmp_ex2.degree(*pit)<<", ";
                
                tmp_ex2 =  tmp_ex2.subs(*pit==1);
          }

      if(r_nu(nu_vec)>r_max) r_max = r_nu(nu_vec);
      if(s_nu(nu_vec)>s_max) s_max = s_nu(nu_vec);
      ID_set.insert(ID(nu_vec));
  MyIntegrals<<"INT["<<t_nu(nu_vec)<<","<<r_nu(nu_vec)<<","<<s_nu(nu_vec)<<","<<ID(nu_vec)<<","<<nu_str(nu_vec)<<"]; "<<endl; 
//* ( "<<tmp_ex2<<" ) "<<endl;         !!!!!!!!!!!!!!!!!!!!!!IMPORTANT out expr!!!!!!
  //      MyIntegrals<<"}];"<<endl;
        cout<< "ID: "<<ID(nu_vec)<<" r "<<r_nu(nu_vec)<<" s "<<s_nu(nu_vec)<<" t "<<t_nu(nu_vec)<< endl;
        nu_vec_lst.push_back(nu_vec);
    }
  std::ofstream prepare_in("prepare.in");
  prepare_in<<"Processes 4;"<<endl<<endl<<"ReduceIDRecursive";
  for(std::set<unsigned int>::reverse_iterator it = ID_set.rbegin();it!=ID_set.rend();++it)
    prepare_in<<" "<<*it;
  prepare_in<<";"<<endl<<endl<<"NrofEqperFile 500;\n\nR2 "<<r_max<<";\nS2 "<<s_max<<";\n\n"
            <<"UseIBP  true;\nUseLI   false;\n\nEpsilon ep;\nDimensionRule d = 4+ep;\nSeries 6;\n";
    cout<<"P_map"<<sp_P_subs_map<<"ddfdfd        "<<distrib_int<<"    terms  "<<distrib_int.nops()<<endl;
}
  catch(string& what)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      cout<<"*********************************************"<<endl;
    }
  //  sp_v(cnst::k1_i3);
}

lst props(ex in_ex)
{
  exset prop_set;
  in_ex.find(bubblesii::prop(wild(1),wild(2),wild(3)),prop_set);
  return set2lst( prop_set);  
}

ex get_no_prop_part(ex in_ex)
{
  in_ex = in_ex.collect(props(in_ex),true);
 try
  {
    cout<<"No pro pin_ex " <<in_ex<<endl;
  if(is_a<mul>(in_ex))
    {
      lst unit_lst;
      for(unsigned int i =0;i< props(in_ex).nops();i++)
        unit_lst.append(1);
      return in_ex.subs(props(in_ex),unit_lst);
    }
  else throw std::logic_error("Expression is not a MUL,but SUM with different sets of propagators");
  }
 catch (std::exception& e)
  {
    cout << "exception caught: " << e.what() << endl;
  }

 
}

ex get_prop_part(ex in_ex)
{
}

void reduze_ng(Graph& G)
{
try
  {
  //  cout<<"reduze expr to work with"<<ex_from_graph(G).normal().expand(expand_options::expand_indexed).simplify_indexed(reduze_sp())<<"  sp  "<<sp_v(cnst::k1_k3).first<<" "<<sp_v(cnst::k1_k3).second<<endl;
    ex in_norm_expr = ex_from_graph(G).simplify_indexed(reduze_sp());
    cout<< "/integral begin:  "<<in_norm_expr<<endl<<endl;
    exset prop_set;
    in_norm_expr.find(bubblesii::prop(wild(1),wild(2),wild(3)),prop_set);
    cout<<in_norm_expr.expand().collect(props(in_norm_expr),true)<<endl<<endl;
    cout<<prop_set<<endl;
  int i_cntr = 0;
  for(const_iterator it = in_norm_expr.begin();it!=in_norm_expr.end();++it)
    {
      i_cntr++;
      cout<< "Integral "<<i_cntr<<"  = "<<*it<<endl;
    }

  // Topo.in filling
  std::ofstream topo_in("Topo.in_");
  topo_in<<"LoopMomenta";
  for(lst::const_iterator it = cnst::km.begin();it!=cnst::km.end();++it)
    if(in_norm_expr.has(*it))topo_in<<" "<<*it;
  topo_in<<";"<<endl<<endl;
  topo_in<<"ExtMomenta";
  for(lst::const_iterator it = cnst::im.begin();it!=cnst::im.end();++it)
    if(in_norm_expr.has(*it))topo_in<<" "<<*it;
  topo_in<<";"<<endl<<endl
         <<"Dimension d;"<<endl<<endl;

  std::stringstream symbols_str,propagators_str;
  symbols_str<<"Symbol";
  ex  N = in_norm_expr.numer().expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  ex D =  in_norm_expr.denom();//.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  cout<<"reduze num  "<<N<<"   denom  "<<D<<endl;  
  Graph::edge_iterator ei,ei_end;
  exmap repls;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    if(D.match(wild(2)*(indexed(G[*ei].gmomentum,varidx(wild(1),4))*indexed(G[*ei].gmomentum,varidx(wild(1),4,true)) -wild(3)),repls )) cout<< "has momentum   "<<G[*ei].gmomentum<<"   "<<repls<<endl;
    else  if(D.match(wild(2)*(indexed(wild(6),varidx(wild(1),4))*indexed(wild(6),varidx(wild(1),4).toggle_variance()) +wild(3)),repls )) cout<< "has (-) momentum   "<<G[*ei].gmomentum<<"   "<<repls<<endl;
    else cout<< "no momentum   "<<G[*ei].gmomentum<<endl;
  symbol P0("P0");
  symbol x0("x0");
  lst P_lst;
  //P_lst.append(P0);
  exmap P_map;
  //P_map[P0] = 1;
  lst x_lst;
  x_lst.append(x0);
  ex factor = 1;
  int cntr = 0;
  string str;
  ex  PX;   // Denominator in form P_i*X_i, where P_i = (q_i^2-m_i^2)
  PX = x0;
  ex PX_symb = x0; // where P_i are symbols
  ex D_symb = 1;
  exset mass_set; //set for symbols storage, no duplicated symbols
  for(const_iterator it = D.begin();it!=D.end();++it)
    {
      if( is_a<GiNaC::numeric>(*it))
        factor *=(*it);
      else
        {
          ex prop_ex;
          exmap prop_repls;
          ex found_mass(wild(1));
          if(it->match(wild()-pow(wild(1),2),prop_repls))
            {
              found_mass = found_mass.subs(prop_repls);
              cout<<prop_repls<<" momentum "<<*it<<endl;
              prop_ex = *it;
              cntr++;
            }
          else if(it->match(wild()+pow(wild(1),2),prop_repls))
            {
              found_mass = found_mass.subs(prop_repls);
              cout<<prop_repls<<" (-) momentum "<<*it<<endl;
              factor*=(-1);
              prop_ex = (*it)*(-1);
              cntr++;
            }
          else throw  string("Bad devision in rational expression");
          //          cout<<prop_ex<<"  wildcard "<<indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true))-pow(wild(2),2)<<endl;
          // extracting mass and momentum from propagators
          mass_set.insert(found_mass);
          exset found_set;
          exmap find_match;
          if(prop_ex.find(indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true)),found_set))
            {
            if(found_set.size() == 1)
              {
                if( found_set.begin()->match(indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true)),find_match))
                  propagators_str<<"Propagator  "<<ex(wild()).subs(find_match)<<" "<<found_mass<<";"<<endl;
                  //  cout<<"momrepls "<<*(found_set.begin())<<"  Propagator  "<<ex(wild()).subs(find_match)<<" "<<found_mass<<";"<<endl;
              }
            else throw string("Bad denominator expression");
            }
          str = "P" + lexical_cast<string>(cntr);
          symbol p_sym(str);
          P_lst.append(p_sym);
          P_map[p_sym] = prop_ex.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
          str = "x" + lexical_cast<string>(cntr);
          symbol x_sym(str);
          x_lst.append(x_sym);
          PX+=x_sym*P_map[p_sym];
          PX_symb+=x_sym*p_sym;
          D_symb*=p_sym;

        }
    }
  topo_in<<"Symbol";
  for(exset::iterator it = mass_set.begin();it!=mass_set.end();++it)
    topo_in<<" "<<*it;
  
  topo_in<<";"<<endl<<endl
         <<propagators_str.str()<<endl;
  N *=factor;

  lst sp_lst(get_sp_lst());
  lst N_has_sp_lst;
  lst PX_has_sp_lst;
  lst irreducible_sp_lst;
  exmap sp_P_subs_map;
  for(lst::const_iterator it = sp_lst.begin();it!=sp_lst.end();++it)
    {
      if(N.has(*it))N_has_sp_lst.append(*it);
      if(PX.has(*it))PX_has_sp_lst.append(*it);
      if(N.has(*it)&&(!PX.has(*it)))irreducible_sp_lst.append(*it);
    }
  for(lst::const_iterator nit = N_has_sp_lst.begin();nit!=N_has_sp_lst.end();++nit)
    {
      // ex (*nit) = N_has_sp_lst.op(0);
  ex PX_sub_part = 0;
  lst eqns;
  for(lst::const_iterator it = PX_has_sp_lst.begin();it!=PX_has_sp_lst.end();++it)
    {
      if((*nit) == *it) 
        {
          eqns.append(PX.coeff(*it,1)==1);
          PX_sub_part += PX.coeff(*it,1)*(*it);           
        }
      else
        {
          eqns.append(PX.coeff(*it,1)==0);
          PX_sub_part += PX.coeff(*it,1)*(*it);           
        }
    }
  eqns.append((PX-PX_sub_part).collect(sp_lst) ==0);


  //  cout<<"P_map "<<P_map<<"    PX      "<<PX<<"   sp__has  "<<irreducible_sp_lst<<"  scalar part  "<<(PX-PX_sub_part).collect(sp_lst)<<"  eqns   "<<eqns<< "    result    "<<lsolve(eqns,x_lst)<<endl;                    
  sp_P_subs_map[*nit] = PX_symb.subs(lsolve(eqns,x_lst));
    }
  
  ex the_int = N.subs(sp_P_subs_map);
  for(lst::const_iterator it = P_lst.begin();it!=P_lst.end();++it)
    the_int/=(*it);
  ex distrib_int = the_int.expand().collect(P_lst,true);


  std::list<std::vector<int> > nu_vec_lst;
  std::set<unsigned int> ID_set;
  unsigned int r_max = 0;
  unsigned int s_max = 0;
  // output file for integrals
  std::ofstream MyIntegrals("afdc_Integrals");
  for(const_iterator it = distrib_int.begin();it!=distrib_int.end();++it)
    {
      std::vector<int> nu_vec;
      ex tmp_ex2 = *it;
      for(lst::const_iterator pit = P_lst.begin();pit!=P_lst.end();++pit)
          {
                nu_vec.push_back(-tmp_ex2.degree(*pit));
                //cout<<-tmp_ex2.degree(*pit);
                //MyIntegrals<<-tmp_ex2.degree(*pit);
                //               cout<<-tmp_ex2.degree(*pit)<<", ";
                // MyIntegrals<<-tmp_ex2.degree(*pit)<<", ";
                
                tmp_ex2 =  tmp_ex2.subs(*pit==1);
          }

      if(r_nu(nu_vec)>r_max) r_max = r_nu(nu_vec);
      if(s_nu(nu_vec)>s_max) s_max = s_nu(nu_vec);
      ID_set.insert(ID(nu_vec));
  MyIntegrals<<"INT["<<t_nu(nu_vec)<<","<<r_nu(nu_vec)<<","<<s_nu(nu_vec)<<","<<ID(nu_vec)<<","<<nu_str(nu_vec)<<"]; "<<endl; 
//* ( "<<tmp_ex2<<" ) "<<endl;         !!!!!!!!!!!!!!!!!!!!!!IMPORTANT out expr!!!!!!
  //      MyIntegrals<<"}];"<<endl;
        cout<< "ID: "<<ID(nu_vec)<<" r "<<r_nu(nu_vec)<<" s "<<s_nu(nu_vec)<<" t "<<t_nu(nu_vec)<< endl;
        nu_vec_lst.push_back(nu_vec);
    }
  std::ofstream prepare_in("prepare.in");
  prepare_in<<"Processes 4;"<<endl<<endl<<"ReduceIDRecursive";
  for(std::set<unsigned int>::reverse_iterator it = ID_set.rbegin();it!=ID_set.rend();++it)
    prepare_in<<" "<<*it;
  prepare_in<<";"<<endl<<endl<<"NrofEqperFile 500;\n\nR2 "<<r_max<<";\nS2 "<<s_max<<";\n\n"
            <<"UseIBP  true;\nUseLI   false;\n\nEpsilon ep;\nDimensionRule d = 4+ep;\nSeries 6;\n";
    cout<<"P_map"<<sp_P_subs_map<<"ddfdfd        "<<distrib_int<<"    terms  "<<distrib_int.nops()<<endl;
}
  catch(string& what)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      cout<<"*********************************************"<<endl;
    }
  //  sp_v(cnst::k1_i3);
}




void bubble(Graph& G)
{
try
  {
  //  cout<<"reduze expr to work with"<<ex_from_graph(G).normal().expand(expand_options::expand_indexed).simplify_indexed(reduze_sp())<<"  sp  "<<sp_v(cnst::k1_k3).first<<" "<<sp_v(cnst::k1_k3).second<<endl;
  ex in_norm_expr = ex_from_graph(G).normal();

  // Topo.in filling
  std::ofstream topo_in("tmp1");
  topo_in<<"LoopMomenta";
  for(lst::const_iterator it = cnst::km.begin();it!=cnst::km.end();++it)
    if(in_norm_expr.has(*it))topo_in<<" "<<*it;
  topo_in<<";"<<endl<<endl;
  topo_in<<"ExtMomenta";
  for(lst::const_iterator it = cnst::im.begin();it!=cnst::im.end();++it)
    if(in_norm_expr.has(*it))topo_in<<" "<<*it;
  topo_in<<";"<<endl<<endl
         <<"Dimension d;"<<endl<<endl;

  std::stringstream symbols_str,propagators_str;
  symbols_str<<"Symbol";
  ex  N = in_norm_expr.numer().expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  ex D =  in_norm_expr.denom();//.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
  cout<<"reduze num  "<<N<<"   denom  "<<D<<endl;  
  Graph::edge_iterator ei,ei_end;
  exmap repls;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    if(D.match(wild(2)*(indexed(G[*ei].gmomentum,varidx(wild(1),4))*indexed(G[*ei].gmomentum,varidx(wild(1),4,true)) -wild(3)),repls )) cout<< "has momentum   "<<G[*ei].gmomentum<<"   "<<repls<<endl;
    else  if(D.match(wild(2)*(indexed(wild(6),varidx(wild(1),4))*indexed(wild(6),varidx(wild(1),4).toggle_variance()) +wild(3)),repls )) cout<< "has (-) momentum   "<<G[*ei].gmomentum<<"   "<<repls<<endl;
    else cout<< "no momentum   "<<G[*ei].gmomentum<<endl;
  symbol P0("P0");
  symbol x0("x0");
  lst P_lst;
  //P_lst.append(P0);
  exmap P_map;
  //P_map[P0] = 1;
  lst x_lst;
  x_lst.append(x0);
  ex factor = 1;
  int cntr = 0;
  string str;
  ex  PX;   // Denominator in form P_i*X_i, where P_i = (q_i^2-m_i^2)
  PX = x0;
  ex PX_symb = x0; // where P_i are symbols
  ex D_symb = 1;
  exset mass_set; //set for symbols storage, no duplicated symbols
  for(const_iterator it = D.begin();it!=D.end();++it)
    {
      if( is_a<GiNaC::numeric>(*it))
        factor *=(*it);
      else
        {
          ex prop_ex;
          exmap prop_repls;
          ex found_mass(wild(1));
          if(it->match(wild()-pow(wild(1),2),prop_repls))
            {
              found_mass = found_mass.subs(prop_repls);
              cout<<prop_repls<<" momentum "<<*it<<endl;
              prop_ex = *it;
              cntr++;
            }
          else if(it->match(wild()+pow(wild(1),2),prop_repls))
            {
              found_mass = found_mass.subs(prop_repls);
              cout<<prop_repls<<" (-) momentum "<<*it<<endl;
              factor*=(-1);
              prop_ex = (*it)*(-1);
              cntr++;
            }
          else throw  string("Bad devision in rational expression");
          //          cout<<prop_ex<<"  wildcard "<<indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true))-pow(wild(2),2)<<endl;
          // extracting mass and momentum from propagators
          mass_set.insert(found_mass);
          exset found_set;
          exmap find_match;
          if(prop_ex.find(indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true)),found_set))
            {
            if(found_set.size() == 1)
              {
                if( found_set.begin()->match(indexed(wild(),varidx(wild(1),4))*indexed(wild(),varidx(wild(1),4,true)),find_match))
                  propagators_str<<"Propagator  "<<ex(wild()).subs(find_match)<<" "<<found_mass<<";"<<endl;
                  //  cout<<"momrepls "<<*(found_set.begin())<<"  Propagator  "<<ex(wild()).subs(find_match)<<" "<<found_mass<<";"<<endl;
              }
            else throw string("Bad denominator expression");
            }
          str = "P" + lexical_cast<string>(cntr);
          symbol p_sym(str);
          P_lst.append(p_sym);
          P_map[p_sym] = prop_ex.expand(expand_options::expand_indexed).simplify_indexed(reduze_sp());
          str = "x" + lexical_cast<string>(cntr);
          symbol x_sym(str);
          x_lst.append(x_sym);
          PX+=x_sym*P_map[p_sym];
          PX_symb+=x_sym*p_sym;
          D_symb*=p_sym;

        }
    }



  topo_in<<"Symbol";
  for(exset::iterator it = mass_set.begin();it!=mass_set.end();++it)
    topo_in<<" "<<*it;
  
  topo_in<<";"<<endl<<endl
         <<propagators_str.str()<<endl;
  N *=factor;

  lst sp_lst(get_sp_lst());
  lst N_has_sp_lst;
  lst PX_has_sp_lst;
  lst irreducible_sp_lst;
  exmap sp_P_subs_map;
  for(lst::const_iterator it = sp_lst.begin();it!=sp_lst.end();++it)
    {
      if(N.has(*it))N_has_sp_lst.append(*it);
      if(PX.has(*it))PX_has_sp_lst.append(*it);
      if(N.has(*it)&&(!PX.has(*it)))irreducible_sp_lst.append(*it);
    }
  for(lst::const_iterator nit = N_has_sp_lst.begin();nit!=N_has_sp_lst.end();++nit)
    {
      // ex (*nit) = N_has_sp_lst.op(0);
  ex PX_sub_part = 0;
  lst eqns;
  for(lst::const_iterator it = PX_has_sp_lst.begin();it!=PX_has_sp_lst.end();++it)
    {
      if((*nit) == *it) 
        {
          eqns.append(PX.coeff(*it,1)==1);
          PX_sub_part += PX.coeff(*it,1)*(*it);           
        }
      else
        {
          eqns.append(PX.coeff(*it,1)==0);
          PX_sub_part += PX.coeff(*it,1)*(*it);           
        }
    }
  eqns.append((PX-PX_sub_part).collect(sp_lst) ==0);


  //  cout<<"P_map "<<P_map<<"    PX      "<<PX<<"   sp__has  "<<irreducible_sp_lst<<"  scalar part  "<<(PX-PX_sub_part).collect(sp_lst)<<"  eqns   "<<eqns<< "    result    "<<lsolve(eqns,x_lst)<<endl;                    
  sp_P_subs_map[*nit] = PX_symb.subs(lsolve(eqns,x_lst));
    }
  
  ex the_int = N.subs(sp_P_subs_map);
  for(lst::const_iterator it = P_lst.begin();it!=P_lst.end();++it)
    the_int/=(*it);
  ex distrib_int = the_int.expand().collect(P_lst,true);


  std::list<std::vector<int> > nu_vec_lst;
  std::set<unsigned int> ID_set;
  unsigned int r_max = 0;
  unsigned int s_max = 0;
  // output file for integrals
  std::ofstream MyIntegrals("tmp");
  for(const_iterator it = distrib_int.begin();it!=distrib_int.end();++it)
    {
      std::vector<int> nu_vec;
      ex tmp_ex2 = *it;
      for(lst::const_iterator pit = P_lst.begin();pit!=P_lst.end();++pit)
          {
                nu_vec.push_back(-tmp_ex2.degree(*pit));
                //cout<<-tmp_ex2.degree(*pit);
                //MyIntegrals<<-tmp_ex2.degree(*pit);
                //               cout<<-tmp_ex2.degree(*pit)<<", ";
                // MyIntegrals<<-tmp_ex2.degree(*pit)<<", ";
                
                tmp_ex2 =  tmp_ex2.subs(*pit==1);
          }

      if(r_nu(nu_vec)>r_max) r_max = r_nu(nu_vec);
      if(s_nu(nu_vec)>s_max) s_max = s_nu(nu_vec);
      ID_set.insert(ID(nu_vec));
  MyIntegrals<<"INT["<<t_nu(nu_vec)<<","<<r_nu(nu_vec)<<","<<s_nu(nu_vec)<<","<<ID(nu_vec)<<","<<nu_str(nu_vec)<<"]; "<<endl; 
//* ( "<<tmp_ex2<<" ) "<<endl;         !!!!!!!!!!!!!!!!!!!!!!IMPORTANT out expr!!!!!!
  //      MyIntegrals<<"}];"<<endl;
        cout<< "ID: "<<ID(nu_vec)<<" r "<<r_nu(nu_vec)<<" s "<<s_nu(nu_vec)<<" t "<<t_nu(nu_vec)<< endl;
        nu_vec_lst.push_back(nu_vec);
    }
  std::ofstream prepare_in("tmp3");
  prepare_in<<"Processes 4;"<<endl<<endl<<"ReduceIDRecursive";
  for(std::set<unsigned int>::reverse_iterator it = ID_set.rbegin();it!=ID_set.rend();++it)
    prepare_in<<" "<<*it;
  prepare_in<<";"<<endl<<endl<<"NrofEqperFile 500;\n\nR2 "<<r_max<<";\nS2 "<<s_max<<";\n\n"
            <<"UseIBP  true;\nUseLI   false;\n\nEpsilon ep;\nDimensionRule d = 4+ep;\nSeries 6;\n";
    cout<<"P_map"<<sp_P_subs_map<<"ddfdfd        "<<distrib_int<<"    terms  "<<distrib_int.nops()<<endl;
}
  catch(string& what)
    {
      cout<<"*********************************************"<<endl;
      cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      cout<<"*********************************************"<<endl;
    }
}




